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We study in this work steady laminar flows in a low density granular gas modelled 
Q ' as a system of identical smooth hard spheres that collide inelastically. The system is 

"^ I excited by shear and temperature sources at the boundaries, which consist of two infi- 

nite parallel walls. Thus, the geometry of the system is the same that yields the planar 
r^ , Fourier and Couette flows in standard gases. We show that it is possible to describe the 

^ • steady granular flows in this system, even at large inelasticities, by means of a (non- 

'"^ I Newtonian) hydrodynamic approach. All five types of Couette-Fouricr granular flows 

M ' are systematically described, identifying the different types of hydrodynamic profiles. 

Excellent agreement is found between our classification of flows and simulation results. 
Also, we obtain the corresponding non-linear transport coefficients by following three in- 
dependent and complementary methods: (1) an analytical solution obtained from Grad's 
13-moment method applied to the inelastic Boltzmann equation, (2) a numerical solution 
of the inelastic Boltzmann equation obtained by means of the direct simulation Monte 
\^ • Carlo method and (3) event-driven molecular dynamics simulations. We find that, while 

PsJ I Grad's theory does not describe quantitatively well all transport coefficients, the three 

^D ■ procedures yield the same general classification of planar Couette-Fourier flows for the 

^f^ ' granular gas. 

o' 

(N 

^^ ' 

1. Introduction 

There have been in the recent years a large number of studies on the dynamics of 

/\^ , granular gases, where 'granular gas' is a term used to refer to a low density system of 

many mesoscopic particles that collide inelastically in pairs. Due to inelasticity in the 

collisions, the granular gas particles tend to collapse to a r est sta te, unless there is some 

kind of energy input. In particular. iGoldhirsch fc Zanettil (|l993l) showed that clustering 



instabilities spontaneously appear in a freely evolving granular gas. Nevertheless, most 
situations of practical interest involve an energy input to compensate for the energy 
loss and sustain, in some cases, the 'gas' condition of the granular system. This type 
of problem has been exte nsively studied, giving r i se to a subfield of granular dynam- 
ics: 'rapid granular flows' ( Jenkins fc Savage! llQSSt IWang. Jackson fc SundaresanI Il996l : 



Goldhirschll2003t lAranson fc Tsimring||2006l) . Furthermore, it has been shown that rapid 



granular flows can attain steady states, some of which, under appropriate circumstances 
and for simple geometries, can give rise to laminar flows, in the same way as a regular 



gas does (see, for instance, the work bv lTii. Tahiri. Montanero. Garzo. Santos fc Duftv 
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2OOIL on Couette granular flows). The question arising ( Goldhirscbll2003h is, what is the 



appropriate theoretical approach to study these granular flows? 

Let us start with classical non-eq uilibrium statistical mechan ics for an ideal gas de- 
scribed by the Boltzmann equation ( Chapman fc Cowling! 1 19701) . As is well known, the 
equilibrium velocity distribution fun ction flr,v,t ) for an ordinary (i.e., elastic) gas is 
the Maxwell-Boltzmann distribution ( Huanelll987l) . For non-equilibrium states, however, 
the solution of the Boltzmann equation is generally not known. On the other hand, in 
some cases, there exist special solutions where all the space and time dependence of 
f{r,v,t) occurs only through a functional dependence on the average fields n (density), 
u (flow velocity) and T ( temperature) associated wi th the conserved quantities (mass, 
momentum and energy) ( Chapman fc Cowling 19701). This typ e of solution is called a 
normal solution of the Boltzmann equation ([Cercignanil 119881 ) . As a consequence, the 
momentum and heat fluxes are also functionals of the hydrodynamic fields and thus the 
balance equations become a closed set of equations for those fields. Ther efore, the normal 
solutions of the Boltzmann equation yield a hydrodynamic description (lHafnll983l ). since 
the closed set of equations is actually for mally similar to the traditional fluid mechan- 
ics equations ( Chapman fc Cowlinall970r ). In practice, what we have got is a transition 
from a microscopic description (based o n the distribu tion function) to a macroscopic 
description (based on the average flelds) (|Hilbertill912l) . 

When the strength of the hydrodynamic gradients is small, the above functional depen- 
dence of the non-uniform distribu tion function / on n,ii,T ca n be constructed by means 
of the Chapman-Enskog method ( Chapman fc CowlinEj|l970l ). whereby / is expressed as 
a series in a formal parameter e: 



f(i) 



/ = .r" + .r>f- + .r>e' + / 



(2),2 



P(3)g3 



(1.1) 



The parameter e indicates the order in the spatial gradients of the average flelds, scaled 
with the inverse of a typical microscopic length unit (mean free path, for instance). If 
terms up to only first order in the gradients are considered (/ ~ /^°^ -I- /^^^e), the mass, 
momentum and energy b alance equations are the wel l known Nayier-S tokes (NS) equa- 
tions of fluid mechanics (IChapman fc Cowlina ll97Cll : ICercignanil 119881 ). This approach 
is accurate for problems where the spatial gradients are sufficiently small. For not so 
small gradients, terms up to secon d order in the gradients need to be considered, and 



we ob tain t he Burnett equations (|Burnettl 119351), used for instance in rarefied gases 



{ Mont ancro . Lopez de Haro. Garzo fc Santosll998tlMontanero. Lopez de Haro. Santos fc Garzd 



1999t lAgarwal. Yun fc Balakrishnanll200ll ). For both NS and Burnett equations, the ex- 
pressions for the fiuxes include a set of parameters called 'hydrodynamic transport coef- 
ficients'. 

Regarding the granular gas, and from a theoretical point of view, it makes sense in prin- 
ciple, due to the system's low density, to derive the dynamics from a closed kinetic equa- 
tion fo£_the_dlstribution function of a single particle, in an analogous way to the standard 
gas (|Goldhirschll2003[ ) : i.e., it is assumed that pre-coUisional velocities are not statistically 
correlated (or, at least, that their correlations are not important). Thus, the correspond- 
ing kinetic equation is analogous to the Boltzmann equat ion but with the modifica- 



tion that inelasticity iii troduces in the collision integral part (jBrev. Duftv. Kim fc Santos 



1998t lGoldhirschll2003 ). We ma y call this modifi e d version of the Boltzmann equation 



'inelastic Boltzmann equation' (iBrev et al 



19981: iGoldhirschI 120031 ). In addition, if we 
assume the existence of a normal solution to the inelastic Boltzmann equation, a hy- 
drodynamic description analogous to that described above for an elastic gas results for 
a granular gas; i.e., transport coefficients and a set of hydrodynamic equations may be 
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derived. This is obviously a question of much interest in the description of transport 
properties of large sets of grains at low density. 

However, due to the coupling between spatial gradient s and inelasticity in steady states 
( Sela fc Goldhirschlll998t ISantos. Garzo fc Duftvl 12004 ). the collisional cooling sets the 
strength of the spatial gradients and thus scale separation might n ot occur (i.e., gradients 
migh t not be small), except in the limit of quasi-elastic collisions (jVega Reves fc Urbach 
20091 ). Therefore, NS or Burnett hydrodynamics would only be expected to work well 
for steady granular flows in the quasi-elastic limit. Nevertheless, some recent works have 
found that a non-Newtonian hydrodynamic description of planar laminar flows, beyond 
Burne tt order, is s till possible for m oderately large spatial gradient s, even for large inelas- 
ticity ( Tii. Tahiri. Montancro. Garzo. Santos fc Duftv;!2001: Santos. Garzo fc Veea Reved 



ty llin. laniri. Montancro. ijarzo. pantos &; JJuitvi/UUi: aantos. uarzo dz Vega neve 
9t lVega Reyes. Santos fc Garz6ll20ld : IVega Reves. Garzo fc Santodboilol ). Actually, 



2008 



it is not surprising that a generalized hydrodynamic description of the Boltzmann in- 
elastic equation works in rapid granular flows, even for moderately large gradients , since 
this is also possible w hen strong gradients occur in elastic gases ( Agarwal et adl200ll : 
Garzo fc Santog|2003[ ). Wc have pointed out previously that this implies that hydrody- 
namics for granular gases is a generalization of classic hydrodynamics for elastic gases. 
Furthermore, a special class of flows has been re cently found in a unified hydrod ynamic 
description valid for elastic and inelastic gases ( Vega Reves et a/.ll2010l uOllof ). Thus, 
the only formal difference between transport theory for granular and ordinary gases 
would emerge not from the limitations due to scale separation but from the possible 
influence of statistical correlations arising from memory effects due to inelasticity. In 
fact, there is a number of works showi ng velocity correlations in systems of inelastic 
particles (for instance, see the work b y McNamara fc Ludind 1998 ;_ Soto fc Mareschall 
200ll:ISoto. Piaseck i fc Marcschal"200l':'Pagonabarr aga. Trizac. van Noiic fc Ernstll2002l : 
Prevost. Egolf fc Urbach 2002: Brilliantov, Poschel. Kranz fc Zippelius 20071 ) and elastic 



particles ([Schlamp fc HathornI 120071 ). This statistical effect would have its origin at the 
more fundamental level of the kinetic equation (the inelastic Boltzmann equation). Put 
in other words, if the Boltzmann inelastic equation is to be valid, hydrodynamic solutions 
for steady gran ular flows arising from it should work, as has been previously shown by dif- 
ferent authors (JAlam fc Nottlll998HTii et a/.|[20oHVega Reves et al\\20m . As a matter 
of fact, the inelastic Boltzmann equation has been used, with good resul ts, as the starting 
point in an overwhelming number of studies on rapid granular flows ( GoldhirschI 12003 ; 
Aranson fc Tsimrindl2006l ). Additionally, good agreement has also been shown, for a va- 
riety of rapid granular flows, between hydrodynamic theory (stemming from the inelastic 
Boltzmann equation) and molecular dynamics results (i n which the velocity statistical 
correlatio ns would be inherently pr esent, see the works bvlPrevost. Egolf fc U rbach 2002: 
Lutsko. Brcv fc Duftv 2002: Dahl. Hrenya. Garzo fc Duftvl l2002t lAlam"fcL uding 20033; 
Montancro. Garzo. Alam fc Ludingl 120061 ). Furthermore, in the case of the special class 
mentioned before, the agr eement of molecul ar dyna mics results with (Grad's) hydrody- 
namic theory is excellent ( Vega Reves et aLl l2010(. 2 011a[ ). 

A considerable amount of work has been devoted to systematic calculations of hydrody- 
namic transport coefficients for granular gas systems, with different d egrees of a pproach 
in th e perturbative solution of the non- uniform distribution function ( Sela fc Gol dhirschI 
1998) : iBrev e^ al\ Il998l : iGoldhirschI l2003t iNott et al\ Il999t I Alam et al\\200^ . However, 
the derivation of non-Newtonian transport coefficients in simple laminar flows has been 
probably not as systematic as for the case of NS transport coefficients. 

The main goal of this paper is the systematic derivation, by means of a non-Newtonian 
hydrodynamic approach, of the steady profiles for laminar granular flows in the simple 
geometry of two infinite parallel walls containing the gas. More specifically, shear and 
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Figure 1. Schematic view of the system subject of study. The granular gas is heated and sheared 
from two infinite parallel walls. Walls are located at j/ = ±/i/2 and have temperatures T± and 
velocities U±, respectively. 



energy are input from the walls (see figure [1} . In the theoretical approach we assume 
that (i) the hydrostatic pressure?) is constant, (ii) the reduced shear rate a (i.e., the ratio 
between the local shear rate and the local collision frequency) is also constant, (iii) the 
shear stress is independent of the granular temperature gradient dyT , whereas (iv) the 
heat flux qy is proportional to dyT. As wc will sec, the resulting classification of profiles 
is formally analogous to the one found for NS hydrodynamics in the quasiclastic limit 
( Vega Reyes fc Urbachll2009l ). except that the constitutive relations arc non- linear. This 
classification is done based on the signs of dy{T^^^dyT) and dyT. As wc will show, both 
signs remain constant throughout the system and are related to the competition between 
viscous heating and inelastic cooling. Moreover, the sign of d^T is also governed by the 
wall temperature difference. In the case of elastic collisions, o nly the viscous heatin g 
effect is present and so dy{T^/^dyT) < 0, which implies d^T < (JGarzo fc Santosll2003l ). 
Therefore, the general classification is only relevant for granular gases and, consequently, 
the case of ordinary gases is embedded as a particular case. 

The hypotheses (i)~(iv) are sensible for a number of reasons. First, they have shown 
a good agreement with computer simulations i n previous work s on Couette granular gas 
flows in the particular case dy{T^^^dyT) < ( Tii et ad 120011 ). I n addition, there exists 
a special class of flows, i ncluding both elastic and inela stic flows (jVega Reves fc UrbachI 
2009HSantos et a/.ll2009l : IVega Reves et aLll20ldEoil(J ). characterized by dy{T^/^dyT) = 



0. This special class defines a surface in the three-parameter space conformed by inelas- 
ticity (represented by the coefficient of normal restitution a), reduced shear rate and 
thermal gradient, as shown in figure [2] It is called 'LTu' surface since this class o f flows 
is characterized by having linear T{ux) profiles ( Vega Reves et all l20ld 12011(3 ). The 
LTu surface splits the parameter space into two regions: the first region (above the LTu 
surface in figure [2] and labelled XTu) corresponds to dy{T^'^dyT) < (i.e., viscous heat- 
ing overcomes inelastic cooling), while the second region (below the LTu surface) has 
dyiT^^'^dyT) > (i.e., inelastic cooling dominates). As we will see, the region below 
the LTu surface can also be split into two sub-regions (labelled CTu/XTy and CTy), 
depending on the sign of dyT, separated by a surface where dyT = 0. The latter surface 
is called here 'LTy' because it corresponds to states where T[y) is a linear function. 
To the best of our knowledge, the regions below the LTu surface have n ot been ex- 
plored before for a 7^ 0, except in the NS description (IVega Reves fc Urbac h 2009). All 
other studies below the LTu surfa ce have been restricted to the plane a = in figure 
2] (see, f or instance, the works by Grossman. Zhou fc Ben-NaimI 119971 : iBrev fc Cubero 
1998; .Brcv. Ruiz-Montero fc Morenoll2000| ). The most prominent result in studies for the 
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Figure 2. Each point of this diagram represents a steady-state Couette-Fourier flow defined 
univocally by the set of parameters ST* (difference between the temperatures at the top and 
bottom fluid layers, divided by the wall separation), a (reduced shear rate) and a (coefficient of 
restitution), the two first ones being determined from the boundary conditions. The surface with 
the label LTu defines the class of states where the temperature T is a linear function of the flow 
velocity Ux, while the surface labelled as LTy (below the LTu surface) deflnes the class with a 
linear profile T{y). Both surfaces intersect in the line representing the uniform shear flow (USF), 
located in the ST* = plane. In addition, the LTu surface contains the line corresponding to 
Fourier flows for ordinary gases (represented by the ST* axis, i.e., a = and a = 1). The point 
ST* — 0, a — and q = 1 (not visible in the diagram) represents the equilibrium state of an 
ordinary gas. Notice that, whereas the LTu surface has points for all values of ST*, the LTy 
surface has an upper bound of ST* which occurs at a = for each a. The LTu and LTy surfaces 
split the space into three regions: XTu, CTu/XTy and CTy (see ij |5.3[) . 



g = plane is perhaps the finding of LTy s tates (IBrev. Cubero. Moreno fc Ruiz-Monterd 
200ll iBrev. Khalil fc Ruiz-Monterol 120091 : iBrev. Khalil fc Duftvll201ll l2012l ). which are 



represented in figure [2] by the intersection curve between the LTy surface and the plane 
a = 0. 

Our purpose is now to extend results obtained in previous works by providing a com- 
prehensive description of granular/clastic Couette-Fourier gas flows, as depicted in figure 
[2] For instance, by determination of t he LTy surface we get to connect the LTy states 
for a = found bv iBrev et all ( 20011 ) with the well known uniform shear flow (USF , 
also referred to as 'simple shear flow', see for instance the works bv ICampbel]||l989[ ). 
within the same theoretical frame. We will follow three complementary routes. First, 
we w ill undertake a theoretical description based on Grad's 13-moment method ( GradI 
19491 ). Second, we will obtain results from two independent simulation methods, the di- 



rect simulation Monte Carlo (DSMC) method, from which a numerical solution of the 
inelastic Boltzmann equation is obtained, and event-driven molecular dynamics (MD) 
simulations, which solve Newton's equations of inelastic hard spheres. As wc will show, 
both simulation techniques support the classification of states mentioned before (and 
sketched in figure [2]). Moreover, the non-Newtonian transport coefficients obtained from 
the approximate Grad solution agree reasonably well with simulations. 

The structure of this work is as follows. In §[2] we describe in more detail the system 
under study and write the corresponding kinetic and average balance equations. For the 
sake of completeness, the solution at the NS level is bricfiy recalled in § [3l Next, the 
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theoretical Grad's solution is derived in § [H In § [5] the assumptions (i)-(iv) referred to 
above are introduced and the associated classification of states is worked out. In § [6] 
we briefly describe the computational methods and compare the simulation results with 
Grad's theory. Finally, we conclude the paper with a summary and discussion in § [T] 



2. Boltzmann kinetic theory and general balance equations 

The system we study is depicted in figure [Ij It is bounded by two infinite parallel 
walls from where we input energy to a granular gas enclosed in between. The energy is 
input by heating (both walls are in general at different temperatures) and, optionally, 
shearing (walls may be moving at different velocities). The granular gas is composed 
by a large number of inelastic smooth hard disks/spheres (inelastic because kinetic en- 
ergy is not conserved during collisions). We consider a set of disks/spheres that is suf- 
ficiently sparse at all times; i.e., the rate at which energy is input is always intense 
enough so that kinetic energy loss in collisions will not cause th e system to 'freeze' or 



'collaps e' (so 'inelastic collapse' does not occur; see for instance iGoldhirsch fc Zanetti 



11993; Kolvin. Livne fc Meersonll2010l) . By sufficiently sparse we mean that we deal with 



a gas in the kinetic theory sense: collisions are only binary and instantaneous (time dur- 
ing collisions is very short compared to typical time between consecutive collisions). We 
consider also that their pre-collision velocities are statistically uncorrelated ('molecular 
chaos' assumption). Therefore, in the absence of external forces, we will assume that 
the veloci ty distribution function of the sys tem obeys the inelastic Boltzmann kinetic 
equation (JBrev et a?.lll998t iBrilhantov fc Poschel 20ol 



^^+v.v\f{r,v;t) = J[v\f,fl (2.1) 

with J being the coUisional integral, whose expression is 

J[v,\f,f]=a''-^ j dv^j d^Q{g.^){g.^)[a-^f{r,v[■,t)f{r,v'^■,t) 

-f{r,vi;t)f{r,V2;t)], (2.2) 

where d is the dimensionality, a is the diameter of a sphere, Q{x) is Heaviside's step 
function, 5 is a unit vector directed along the line joining the centers of the colliding 
pair, g — Vi — V2 is the relative velocity, and {t'i,f2} and {v[,v'2} arc post-coUisional 
and pre-colhsional velocities respectively. As we see in (|2.2p . J[vi\f,f] depends on the 
parameter a, whic h characterizes inelasticity in the c ollisions and is called coefficient of 
normal restitution ( Brev et a/.lll998tlGoldhirschll2003l ). The (restituting) collisional rules 



for a pair of colliding inelastic smooth hard disks/spheres is 

v'2 = v2 + ^{l + a-'){^■g)^. (2.3) 

The first d + 2 velocity moments of f{r,v,t) define the number density n{r,t), the 
flow velocity u{r,t) and the granular temperature T{r,t) as 



dvfiv), (2.4) 

nu = I dvvf{v), (2.5) 
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nT='^JdvV'f{v), (2.6) 

where V = v — u is the peculiar velocity and in is the mass of a particle. 

Mass, niomeiitum and energy balance equations are obtained by multiplying both sides 
of (|2.ip by 1, V, v"^ and integrating over velocity. The results are 

Dtn = -nV-u, (2.7) 

Dtu = —V ■ P, (2.8) 

mn 

DtT + CT = -^{P:Vu + V-q). (2.9) 

an 

In the above equations, £)( = 9t + m • V is the material derivative, 

P = m dvVVfiv) (2.10) 



is the pressure tensor, 

is the heat flux vector and 



q=!^ldvV'Vf{v) 



(2.11) 



^^~Mf /d^^'^t^l/'/] (2-12) 



is the cooling rate characterizing the rate of energy dissipated due to collisions. 

Next, we consider the steady base states that may be generated from energy input in 
our geometry. Independently of the nature of the boundary conditions, and if there is no 
pressure drop source or gravitational field in the horizo ntal directions (which may gener- 
ate Poiseuille flows; see for exa mple the recent works bv lTii fc Santosll2004 ISantos fc Tii 



2006t lAlam fc Chikkadil 120101) , the spatial dependence of these steady base states will 



occur only in the coordinate y, perpendicular to both walls (we call it vertical direction). 
Moreover, the flow velocity is expected to be parallel to the walls, i.e., u{y) = Ux{y)ex. 
Consequently, the Boltzmann equation (|2.1I) for these reference steady states can be 
rewritten as 

' dy 
and the balance equations have the simple forms 

^-0, ^^0, (2.14) 



vl^.^JifJ] (2-13) 



P^y^ + ^]-CT. (2.15) 



dn \ ^^ dy dy 

Due to the symmetry of the problem, all the off-diagonal elements of the pressure tensor 
different from Pxy vanish and, in principle, the two shear-flow plane diagonal elements 
{Pxx and Pyy) are different whereas the remaining d— 2 diagonal elements orthogonal to 
the shear-flow plane are equal. The latter property implies that Pxx + Pyy + {d— 2)Pzz = 
dp, where p = nT = d~^TiP is the hydrostatic pressure. 
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3. Navier Stokes description 

The balance equations (|2.14p and (|2.15p are exact and do not assume any particular 
form, for the constitutive equations. However, they do not constitute a closed set of 
equations for the hydrodynamic fields. 

The simplest approach to close the problem is provided by the NS c onstitutive equa- 
tions, which, in the g eometry of the planar Couette-Fourier flow read (jBrev et aLlll998 : 
Brev fc Cuberoll200ll) 

P.x = Pyy = Pzz = P, (3.1) 

Pxy = -??o'7ns(")-^' (3-2) 

q. = 0, (3.3) 

(3.4) 



^ ^ , .dT T ^ dn 

q.y = -AoKns(")-^ ^o-/^Ns(")^- 



In equations p.2p and p.4p . 

7?o = \/^^c^Ada-('^-i\ 



Ki = 



d + 2 



r{d/2)TT-— 



(3.5) 



is the NS shear viscosity for elastic gases ( Gradl 1 19491 : 1 Chapman fc Cowling|l97Cll ) and 

d{d + 2) CA 7?o 



An = 



(3.6) 



2{d — 1) Cri rn 

is the NS thermal conductivity for elastic gases (|Gracilll949l : IChapman fc Cowlinall97Cl[ ). 
In equations p.Sp and p.6p . the factors c,, and cx take the values c,, = 1.0 22, c\ = 1.029 
for hard disks (d = 2) and £ „ = 1.016, cx = 1.025 for hard spheres (d = 3) (|Burnettill935 



Chapman fc Cowling Il970l ). Finally, 77^3, Kjjg and //Jjg are the reduced NS transport 



coefficients of a dilute granular gas, whose expressions are given in Appendix |^ In 
equations (|Aip - (|A3p . 

C{a)^^{l-a') (3.7) 

represents the ratio between the cooling rate ( and an effective collision frequency defined 
as 

= P. 

~ Vo' 

Note that ly cx nT^/"^ and thus it depends on y. 

Now we combine the NS constitutive equations with the three balance equations (|2.14p 
and (|2.15p . First, the exact property Pyy = const, together with equation p.ip . implies 
that the hydrostatic pressure is uniform. Next, the exact property Pxy = const, together 
with equation p.2p . implies that the product rj^dux/dy = const. These two implications 
can be combined into a — const, where 



(3.8) 



1 dux 
V dy 



(3.9) 



is the reduced shear rate. Finally, we consider the energy balance equation p.l5p . First, 
since p ~ const, equation p.4p can be rewritten as 



qy 



-AoA 



dT 

Ns(")-5r' 



dy' 



-^NS — ^ 



NS 



A^NS- 



(3.10) 



Steady base states for granular hydrodynamics 



Next, using the properties Pxy = const, p = const and a = const in equation ()2.15|) . one 
has v~^dqy/dy = const. This, together with equation p.lOp yields 



1 d fldT\ „ , , 



where 



7ns (a, a) 



1 rj*^^{a)a^ - Ida) 



(3.11) 



(3.12) 



did + 2) a;,s(«) 

Therefore, the NS description, as appHcd to the Couettc-Fourier flow, predicts that the 
hydrostatic pressure p = nT, the reduced shear rate (|3.9p and the second order derivative 
(v~^d„y'T are uniform. A detailed account of this NS description was presented by 
IVega Reves fc UrbachI (|2009l ). 



4. Non-Newtonian description: Grad's 13 moment method 

The results derived in § [3] are restricted to small spatial gradients. Thus, they do not 
capture non-Newtonian effects, such as normal stress differences (i.e., Pxx 7^ Pyy 7^ Pzz) 
and a non-zero component of the heat flux orthogonal to the thermal gradient (i.e., 
qx 7^ 0). Those effects are expected t o be present in the solu tion of the Boltzmann 
equation beyond the quasi-elastic limit (jSela fc GoldhirschlllQQSI ). 



The aim of this section is to unveil those non-Newtonian properties by solving the set of 
moment equ ations derived from the Boltzmann equation by Grad's 13-moment method 
(|Gradlll949l ). In this method, the velocity distribution function / is approximated by the 
form 



/ ^ ./o 1 + 



2nT2 



[P^, - p5.,,)V,V, + 



mV^ 



d + 2 \ 2T 



d + 2 



V q 



where 



/o 



/ ra 



d/2 



-mV'^/2T 



(4.1) 



(4.2) 



\2t:TJ 

is the local equilibrium distribution. The number of moments involved in equation (|4.ip 
is d{d + 5)/2 -I- 1, which becomes 13 in the three-dimensional case. The coefficients in 
Grad's distribution function have been obtained by requiring the pressure tensor and 
heat flux of the trial function (|4.ip to be the same as those of the exact distribution /. 
The Grad distribution (j4.ip can be interpreted as the linearization of the maximum - 
entropy distribution constrained by the first d{d + 5)/2 + 1 moments ( Kremeij 120101 ) . 
From that point of view, it is not guaranteed a priori that it is quantitatively accurate 
for large deviations from the local equilibrium distribution. Moreover , an extra isotropic 
term associated with the fourth velocity moment can also be included (jSela fc Goldhirsch 



19981) . However, here we consider the minimal version of Grad's method, restricting the 
number of non-Maxwellian parameters to the stress tensor and the heat flux vector, since 
extra terms do not significantly increase accuracy. 
According to the approximation (|4.ip . one has 



dvV.VjVkf 



1 



d + 2 



{qiSjk + qjkk + qk5ij); 



(4.3) 



dvV'^V^VJ,f 



p 

nm 



-Pij - pSi 



(4.4) 
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In add ition (JBrev et aLllQQSHBrev fc Cuberol200ll : lGarz6 fc Montanerol2002l : l\^a Reves et d. 



2011aD, 



dv V,VjJ[!, /] ^ -u [/3i {P,j - p5,i) + CP^^ > 



I^jdvV'VJif,/]^- 






(4.5) 



(4.6) 



where, as usual, terms non-linear in P^- — pSij and q have been neglec ted. On the other 
hand , the quadratic term s have been retained in some other works ( Herdegen fc Hesd 
ll982l : lTsao fc Kochlll995[ ). In equations (|4.5p and (|4.6p . the collision frequency v is given 
by p.Sp (and taking into account equation p.Sp ) with c,, = 1. Also, C* = C/^j Pi and /32 
are given by equations p.7p . (|A4p and (|A5p . respectively. 

The relevant moments in our system are p, T, Mj,, P^jy, Pa;^, Pyj/, Qx and g^. The 
exact balance equations (|2.14[) and (|2.15p are recovered by multiplying both sides of 
equation (j2.13p by 14, Vy and V"^ and integrating over velocity. In order to have a closed 
set of differential equations, we need five additional equations, which are obtained by 
multiplying both sides of equation (j2.13p by VxVy, V^, Vy , V'^Vx and V^Vy and applying 
the approximations (|4.3p - (|4.6p . The results are 



d + 2 



:9, 



^yy^s^x 



Wi + C)P. 



(4.7) 



d + 2 



dsQy + 2PxydsUx = -/?! (Pxx ^P)- C P^X , 



d + 2 



dsQy = -Pi [Pyy - p) - CPyy, 



d + 4„ /T 



d + 4 
d + 2 



qydsUx 



d-1 



-P2qx, 



ds 



T 
m 



rf + 4 



Pyy P 



2 Fi d-1 

QxOsUx = -^P2qy, 



d + 2 



where we have introduced the spatial scaled variable s{y) by 

ds = v{y)dy. 



(4.8) 

(4.9) 

(4.10) 

(4.11) 

(4.12) 



Note that dsj ^j2T{y)/m measures the elementary vertical distance dy in units of the 
(nominal) mean free path \j2T{y)/m/i'{y) . Therefore, the scaled variable s{y) has di- 
mensions of speed. Its limit values are deduced from integration of (|4.12p . taking into 
account that the limit values of y are y = +h/2. 

It must be stressed that in equations (|4.7p - (|4.1ip the only assumptions made are the 
stationarity of the system, the geometry and symmetry properties of the planar Couette- 
Fouricr flow and the applicability of Grad's method. 

The exact momentum balance equations p.l4p imply that Pxy — const and Pyy = 
const. Moreover, if one assumes that p = const, equation (j4.9p yields dsQy = const. Next, 
the exact energy balance equation (|2.15p implies that the reduced shear rate a = daUx 
defined by equation (|3.9p is also constant (recall that Q* = Q/v = const). Taking all of 
this into account, we get that dsqx ~ const and Pxx = const from equations (|4.7p and 
(|4.8p . respectively. Finally, equations ()4.10p and ()4.1ip imply that both qx and qy are 



proportional to the thermal gradient dgT. As a consequence, d^T — const. 

Since the pressure p, the shear stress Pxy and the shear rate a = v^^dyUx are constant, 
it follows that the ratio Pxy / ilodyUx is also constant (recall that 770 = p/v)- That ratio 
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defines a (reduced) non-Newtonian shear viscosity coefBcient 77* (a, a) by 

, dux 
dy 



P,y^-TloV*{a,a)^. (4.13) 



Analogously, the fact that qy ex dsT, together with the relationship Aq cx; p/v, allows us 
to define a (reduced) non- Newtonian thermal conductivity coefficient A* (a, a) by 

dT 
qy^-\o\*{a,a) — . (4.14) 

Equations (j4.13p and (|4.14p can be seen as generalizations of Newton's and Fourier's 
law, equations (|3.2p and p.lOp . respectively, in the sense that the reduced transport 
coefficients rj* and A* are non-linear functions of the shear rate a and thus they differ 
from the NS coefficients r/Jjg and AJjg of a granular gas ( Brev et al\\l99S ). It is important 



to not e tnat, due to tne coupling between coili sional cooling and gradients m steady 
states ( Brev fc Cuberdll998t ISantos et al\\2004} . the generalized transport coefficients 



do not reduce to the NS ones in the absence of shearing (a = 0). In fact, at equal wall 
temperatures and in the absence of shearing, an autonomous thermal gradient appears 
in the system that is controlled by inelasticity only, so that A* (a, 0) differs from the NS 
quantity AJjg(a). 

It is interesting to remark that, among the hypotheses (i)-(iv) described in § [TJ only 
the p = const hypothesis is needed in the framework of Grad's set of equations. 

Apart from the generalized coefficients ry* and A*, departures from Newton's and 
Fourier's laws are characterized by normal stress differences and a component of the 
heat flux orthogonal to the thermal gradient. These effects are measured by the (re- 
duced) directional temperatures 

9x{a, a) = — , 9y{a,a)^^, (4.15) 

p p 

and by a cross conductivity coefficient 0* defined as 

dT 

qx = \^cp*{a,a) — . (4.16) 

dy 

Equation (|4.15p is consistent with the fact that the diagonal elements of the pressure 
tensor (i.e., the normal stresses) are uniform, while equation (|4.16p is consistent with 
dsQx = const. The parameters Ox and 6y account for the distinction between the diagonal 
elements {Pxx and Pyy) of the pressure tensor from the hydrostatic pressure p = [Pxx + 
Pyy + {d—2)Pzz]/d. Moreover, cj)* characterizes the presence of a heat flux component qx 
induced by the shearing. These three coefficients are clear consequences of the anisotropy 
of the system created by the shear flow. Note that, by symmetry, the coefficients 77*, A* 
and 9i are even functions of the shear rate a, while (jf is an odd function. 

Inserting equation (|4.13|) into the (exact) energy balance equation (|2.15p . it is straight- 
forward to obtain 

Idqy d{d + 2)^^^ ^ , ^ 

--^=P ^_^ A*(a,a)7(a,a), (4.17) 



with 



( ^_ d-l rj*ia,a)a^-ICH ,, ^., 

'y(a,a) — —y- — -7 r—. ^ . (4.18) 

'^ ' ^ d{d + 2) A*(a,a) ^ ^ 
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Using equation ()4.14|) . equation ()4.17p yields 

-|-f-?-)=-2"^^("'«)- (4.19) 

V oy \v oy J 

The technical steps needed to derive the transport coefficients r/*, A*, Ox, dy and (j)*, 
as well as the thermal curvature parameter 7, in the framework of Grad's method are 
worked out in Appendix |B] 

In summary, we have shown that Grad's 13-moment method to solve the Boltzmann 
equation is consistent with the general assumptions made in § [TJ Moreover, explicit 
expressions for the generalized non-Newtonian transport coefficients are derived. On the 
other hand, given the approximate character of Grad's method, a more quantitative 
agreement with computer simulations is not necessarily expected. 



5. Generalized non-Newtonian hydrodynamics 

5.1. Basic hypotheses 

Sections[3]and|3]show that the exact balance equations (j2.14p and (j2.15p allow for a class 
of base-state solutions characterized by the following features: 

• (i) the hydrostatic pressure p is uniform, 

• (ii) the reduced shear rate defined by equation (|3.9p is uniform, 

• (iii) the shear stress Pxy is a non-linear function of a but is independent of the 
thermal gradient dyT and 

• (iv) the heat flux component qy, properly scaled, is linear in the reduced thermal 
gradient but depends non-linearly on the reduced shear rate a. 

As shown before, in the NS description properties (i)-(iv) are a consequence of the 
constitutive equations themselves, while in the Grad description one only needs to assume 
point (i) and then the other three points are derived. 

It is important to remark that hypotheses (iii) and (iv) are fully consistent with the 
Burnett-order constitutive e quations in the Couette-Fo urier geometry; taking into ac- 
count the general structure ( Chapman fc Cowling||l97Cll ) of the Burnett contribution to 



(2) (2) 

the shear stress, Pxy , and to the heat flux, g^ , it is straightforward to check that 

Pxy = Qy = if WiUj = dyUxSiySjx, VjT = dyTSiy and V^p = dypSiy. 

The aim of this section is to assume the validity of hypotheses (i)-(iv) in the bulk 
domain of the system (i.e., outside the boundary layers) and analyze the different classes 
of base states that are compatible with them. In doing so, we are assuming that the 
Boltzmann equation admits for solutions which, in the bulk domain of the system, are 
essentially in agreement with (i)-(iv), beyond the NS or Grad's app roximations. Pr evious 



results obtained for ordinary (jGarzo fc Santosll2003l ) and granular (JTii et aZ.ll200ir ) gases 
support the above expectation. 

Assumptions (iii) and (iv) can be made more explicit by Eqs. (j4.13p and (|4.14[) . re- 
spectively, where the generalized transport coefficients 77* (a, a) and A* (a, a) have not 
necessarily the explicit forms provided by Grad's solution. The same can be said about 
equations (|4.15p and (14.161) . Moreover, from the energy balance equation (|2.15p one can 
again derive equations (|4.17[) - (|4.19p . provided that the possible spatial dependence of 
the ratio (* = Qjv due to higher-ord er gradients is di scarded. This ass umption is sup- 
ported by kinetic theory calculations ( Brev et aLlll998l ) and simulations ( Tii et a/.ll2001 



Astillero fc Santosll2005[) . 
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According to the assumption p = uT = const, the colhsion frequency defined by 
equation p.Sp has the exphcit form 



^jn-l/2 J^^ 



pa 



d+i 



mCjjAd 



(5.1) 



and thus equation (|4.17p imphes that the product T^^^dyQy is uniform. Moreover, the 
sign of dyQy is determined by that of the coefficient 7. Equivalently, in view of equation 
(j4.19p . the parameter 7 has a direct influence on the curvature of the thermal gradient. 
We see from equation (|4.18p that the main difference between 7 for elastic and inelastic 
gases is the absence or presence of the term proportional to C* , respectively. In both cases 
(i.e., C* = or C* > 0), 7 is constant. On the other hand, while 7 is positive definite in 
the elastic case, its sign results from the competition between viscous heating (?7*a^) and 
inelastic cooling {dC,* /2) in the inelastic case. As a consequence, as we will show below, 
inelasticity s pans a more gener al set of solutions, which includes the clastic profiles as 
special cases (jVega Reves fc Urbach 2009i) . 



5.2. Properties of the hydrodynamic profiles 

In terms of the scaled spatial variable s defined by equation (|4.12l) , equations p.9p and 
(|4.19p take the following forms 

I, (5.2) 



ds 






~2m'y{a, a) 



(5.3) 



From equations (|5.2p and ([5 
terms of the scaled variable: 



it is straightforward to obtain analytical solutions, in 
u.^{s)=as + C, (5.4) 



T{s) = -mj{a,a)s^ + As + B, (5.5) 

where A, B, C are integration constants. Please note that integration of the differential 
equations (|5.2p and (|5.3p is done independently of the nature of the boundary conditions. 
We may set C = by a Galilean transformation. The constants B and A represent the 
values of T and dgT, respectively, at a reference point s = 0. Therefore, since it is always 
possible to choose the point s = within the physical region, henceforth we can take 
B > without loss of generality. Note that equations ()5.4p and ()5.5p imply that T is also 
quadratic when expressed as a function of Ux or, equivalently. 



dui 



-2m 



7(a,a) 



(5.6) 



Taking into account the definition of s and equation (|5.ip (with K = const), we may 



write the derivative 9^T in the natural variable y in terms of dgT and d^T as 



dy'^ ds 



y-1/2 



dT 



K^T- 



T 



9? 



dT 
ds 



By using equation (|5.5p . one gets 



'^^K^T-^Ha,a), 



(5.7) 



(5.8) 
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where $ is also uniform and is defined by 

1 



$(a, a) = -2inB-f{a, a) - ^A^. (5.9) 



'^^--W^^('^ + 2mT7)- (5.12) 



In the same spirit as in equation (|5.6p . the parameter $ can be conveniently expressed 
as 

d^T IfdTV <^{a,a) 

In contrast to 7, the quantity <&, which measures directly the curvature of the thermal 
profile, is determined not only by the shear rate and the inelasticity, but also by the 
temperature boundary conditions through B and A. Similarly, from the identity dyT = 
KT~^/^dsT and equation (j5.5p . it is straightforward to obtain 

/ dT\ ^ 

T(— J = -2K^{'^ + 2mT-f). (5.11) 

This implies that i> is upper bounded: $ ^ ~2mT^. For 7 > 0, one has $ ^ —2mT^axJ, 
while $ ^ 2TOTi„in|7| for 7 < 0. Here, Tmax and T^^ are the maximum and minimum 
values, respectively, of the temperature in the system. Another interesting consequence 
of equation (|5.1ip is that, according to the constitutive equation (|4.14p . qy is a linear 
function of T: 

d^(d + 2) ^p^A*' 
2{d-iy 

The same relationship is obtained for q^, except that A* is replaced by 0*. 

Since both 7 and $ are constant across the system, equations (|5.6p and (|5.8p imply that 
neither T{ux) nor Tiij) exhibit a curvature change, i.e., they do not possess an inflection 
point. On the other hand, this is not necessarily so for the velocity profile Ux{y)- To 
clarify this point, note that, according to equations p.9p and (|5.ip . 

Thus (assuming a > 0), Ux{y) is convex (concave) in the spatial regions where the temper- 
ature increases (decreases). In case the temperature presents a minimum or a maximum 
at a certain point inside the system, the flow velocity presents there an inflection point. 
In the derivation of equation (|5.13p no use of the form of the temperature profile has 
been made. On the other hand, taking derivatives on both sides of equation (|5.13p and 
using equations (|5.8p and (jS.lip . one obtains 

— ^ = -K^aT^'^''^ (2$ + 3mT7) . (5.14) 

Therefore, similarly to T{dT/dy)^ and qf, T'^^^d'^Ux/dy'^ is a linear function of temper- 
ature. 



Eq uations (j5.2p - (|5.14p also apply in the NS hydrodynamic description (jVega Reves fc Urbach 



tiqu 



20091 ) ■ except that ri*{a,a), A* (a, a) and 7(0;, a) are replaced by their NS counterparts 



^Ns('^)' '^Ns('-'^) ^^'^ 7Ns("ja), respectively (see §[3]). While ?/ns('^) ^^'^ '^ns('^) ^^'^ inde- 
pendent of the shear rate, one sees from equation (13.12^ that 7ns (ct, a) is a linear function 
of a^. 
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5.3. General classification of states 

In a previous work (jVega Reves fc UrbachI 120091 ) . the complete set of steady-state solu- 
tions based on the signs of the parameters 7 and $ was described in the framework of NS 
hydrodynamics. It was shown in that work that the analytical expressions of the tem- 
perature and flow velocity profiles depend on the signs of these two parameters. Thus, 
each possible combination of signs of 7 and $ yields a different class of constant pressure 
laminar flows. Now, we can perform the same analysis in the non-Newtonian regime and 
find the same set of classes of steady base states. 
It is convenient to define the following constants 



|$| 



,,2 _ 



l$l 



wTf 



1/2 



to 



_ UJ-LQ 



So 



(5.15) 



2K ' 2m7 

set the natural scales for T, u^ and y. 



277i|7|' 2to^7^ ' 

As we will sec below, the constants Tq, w and 

respectively. According to the signs of 7 and <f>, the following cases are possible: 

(1) 7>0. 
This case [see equation (|4.18|) ] corresponds to states where viscous heating is larger than 
coUisional cooling. Therefore, t his class exists only in the presence of shearing (a ^ 0) 
and inelasticity is not required (JTii et aLll200ll ). Note that, according to equation ()5.9|) . 
7 > implies 

$<0. (5.16) 

From equations ()5.3|) and ()5.6p . T[s) and, equivalently, T{ux) are convex. We will refer 
to this class as XTu. Also, from equations (|5.8p and (|5.16p we conclude that the profile 
T{y) is convex as well. Moreover, equation (j5.12p shows that qf {i = x, y) decreases with 
increasing temperature. 

Making use of the definitions (|5.15p in equation (|5.5p . the quadratic function T[s) can 
be written as 

\ 2' 

s - So 



T{s) = To 



1 



w 



(5.17) 



Since dy = K ^T^/^ds, the relationship between the true and scaled space variables is 



y^yo + ^o 



So 



So 



So 



(5.18) 



Eliminating s between equations (j5.17p and (j5.18p one gets T{y) in implicit form: 

\y-yo\ =4 




sm 



T 
1 

To 



(5.19) 



Equation (j5.18p also provides the velocity profile u^iy) in implicit form just by replacing 
s by u^/a: 



y = yo + ^o 



uo 



uo 



aw 



Uo 



aw 



(5.20) 



where uo = aso ■ A similar replacement in equation (j5.17p yields T as a function of Ux ■ 
In the above equations Sq and yo denote the point where the temperature reaches its 
maximum value T = Tq. This point may be inside the system (i.e., |yo| ^ h/2) or outside 
the system. In the latter case, the maximum corresponds to a continuation of T(y) into 
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the external region |yo| > h/2. The physical condition T(y) > implies the domains 

|s-soKu;, |j/-2/oK|4. (5.21) 

Although the hydrodynamic profiles in terms of the s variable are quite simple [see 
equations ()5.4p and (|5.5|) ]. equations (|5.19|) and (|5.20|) show that the dependence of T 
and Ux on the real space variable y is highly nonlinear. A similar comment applies to the 
cases discussed below (except in the cases LTu and LTy, where the profiles are simpler). 

(2) 7 = 0. 
Now viscous heating exactly equals collisional cooling. As a consequence, T{s) and 
T(u:^) are linear functions. For this reason, we fo rmerly referred to this class as LTu 
( Santos et aLll2009l : IVega Reyes et a/.ll20ld l2011al ). Moreover, the heat flux is uniform 
[see equation ()4.17p ]. 
Two possibilities for $ are found: 

(2.a) $ < 0. 

From equation (|5.9p . A^ = 2|$| 7^ and the profiles are 



T{s) = As + B, 



(5.22) 



Ux{y) = -^ 



T{y) 



2/3 



-AK{y ~ yo) 



-AK{y - yo) 



aB 

~A' 



2/3 



(5.23) 



(5.24) 



Here yo represents the mathematical point where Tiij) — > 0. Obviously, positivity 

of T{y) requires y > yo if A > and y < yo if A < 0. It is possible to prove that 

equation ()5.19|) reduces to equation (|5.24|) in the limit 7—^0. 

Notice that, from equation (j4.18p . 7(0,0) = is fulfilled for a threshold shear rate 

aLTu(a), whose specific value (for a given a) requires the knowledge of ?/* and C*. 

In the special case of elastic collisions (a = 1, i.e., (,* = 0), 7 = implies aj-p^ = 0. 

This corresponds to the conventional Fourier fiow of an ordinary gas. 

(2.b) $ = 0. 

This implies A = 0, so the temperature is uniform and the heat flux vanishes. In 

this case s is a linear function of y and thus equation (j5.4p yields 



Ux{y) = avy 



(5.25) 



with u — const. This stat e is the well-kn own uniform (or simple) shear fiow (USF; 
see, for instance, work bv ICampbelllll989l ) . Note tha t here the USF is not g enerated 
by the usual Lees-Edwards boundary conditions ( Lees fc Edwardslll972h but by 
thermal walls in relative motion. The USF needs again the condition a = aLTu(Q^)- 
Notice that a = \ gives only the trivial equilibrium state of an elastic gas. 
(3) 7 < 0. 
In this wide class, inelastic cooling overco mes viscous heating. Th erefore, collisions must 
be inelastic and shearing is not required ( Brev &: Cuberolll998l ). A negative 7 implies 
a concave curvature of T{s) and T(ux), qf being an increasing (linear) function of T. 
According to equation (j5.9|) . we find now three possibilities for the curvature of the 
temperature profile T{y): 
(3.a) $ < 0. 
In this subclass, henceforth referred to as CTu/XTy, T{y) is a convex function. The 



Steady base states for granular hydrodynamics 



17 



profiles are 



Tis) ^ To 



S" So 



(5.26) 



y^yo + ^o 



So 



s - So 



l-ln 



s - So 




s- So 



- 1 



\y-yo\ = 4 



%{^^^)-m 




(5.27) 
(5.28) 



In equations (|5. 26^ - 1)5. 28p so and yo denote tlie mathematical point where the tem- 
perature reaches its formal minimum value T = —Tq. This point must obviously lie 
outside the system (i.e., \yo\ > h/2). The physical condition T{y) > implies that 



|s-so| ^ w, \y~yo\ ^ -^^o- 



(5.29) 



(3.b) $ = 0. 

This case corresponds to a linear function T{y). Thus, we call this class LTy. From 
equation (|5.9p we have B = A^ /Am\^\ and the profiles are simply 

T{s) =m|7|(s-so)^, 

/ \ 1/2 

so+f^^l [y^yof' 



Ux{y) = a 



vHtI, 



(5.30) 
(5.31) 



T{y) = 2K./^\{y - yo), (5.32) 

where, without loss of generality, we have assumed T{h/2) ^ T{—h/2). Similarly 

to the LTu case, Sq and yo represent the point where T — > 0. Thus, one must have 

y > yo- It is straightforward to reobtain equation (|5.32[) from equation (|5.28|) in the 

limit $ — >■ 0. Note that in the LTy class of states qf /T is constant [see equation 

(I02)) ]. 

If we denote by 

1 AT 
ST*^—^— (5.33) 

the reduced applied gradient, where AT = T{h/2) — T{—h/2), then the LTy flow 
requires a transitional value given by 



sn 



LTy (a, a) = 2^/\"f{a,a)\. (5.34) 

Note that, becau se of expected temperature jumps at the walls (|Lunlll996l : lGalvin et all 
2007tlNottir201l[) . T{±h/2) ^ T±. Moreover, by T{±h/2) here we mean the extrap- 
olation to y = ±ft./2 of the bulk temperature profile, which might differ from the 
respective temperatures of the fluid layers adjacent to the walls, due to boundary- 
layer effects. 

As we will show below, if 7 < 0, I7I always increases with decreasing shear rate a, 
and thus (5T£rp (a, a) has an upper bound at a = given by 



<5rL*Ty(a,a)s:2V|7(a,0)|. 



(5.35) 
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Label 


sign(7) 


sign($) 


Shearing 
needed? 


Inelasticity 
needed? 


ns),T{u,) T{y) 


qlmdiT) 


XTu 
LTu 


+ 



— 


Yes 
Yes* 


No 
Yes* 


Convex Convex 
Linear Convex 


Decreasing 
Constant 


USF (LTu) Vest Yes^ 
CTu/XTy - - No Yes 
LTy - No Yes 
CTy - + No Yes 
* Except for the Fourier flow of an ordinary gas (a = 
^Except for the equilibrium state of an ordinary gas 


Constant Constant 
Concave Convex 
Concave Linear 
Concave Concave 

= 0, Q = l). 

(a = 0, ST* = 0, a = 


Zero 
Increasing 
Increasing 
Increasing 

0- 




Table 


1. Classification of Couette-Fourier flows (see text). 





The LTy state has been studied previously (JBrev et aZ.ll200ll l2009l l201ll l2012l ) in 
the absence of shearing (a = 0). 

In equation ()5.34p it is implicitly assumed that the shear rate a is a free parameter. 
Reciprocally, given an imposed gradient 6T* ^ 2^y\^{a, 0)|, it is always possible to 
find a certain value of the reduced shear rate, auTy{a,ST*), such that 



7(a,aLTy(a,<5T*)) = -l(5T*)^ 



(5.36) 



Since I7I is a decreasing function of a, it is obvious that aLTy increases with decreasing 
ST*. Therefore, the maximum value occurs at ST* = (i.e., 7 = 0), which coincides 
with flLTu (sec figure [2]). In other words, 



auTyio^JT*) i^ aLTu(a)- 



(5.37) 



In fact, the case aLTy = iltu corresponds to the USF state. 

(3.c) $ > 0. 

In this class, T{y) is a concave function and so we call this class CTy. The resulting 

profiles are 

2' 
/ .S — Sn \ 
1 



T(s) = To 



So 



y = yQ + ^Q 



\y-yo\ 




sinh 



So 



sinh"' W— - 1 



T 



(5.38) 
(5.39) 
(5.40) 



where sq and j/o denote the point where the temperature reaches its minimum value 

T = To. 

The main features of the six classes of flows described above arc summarized in table [Ij 

Note that these six profile types have been obtained independently of the specific details 

of the boundary conditions. Once they are specified (and they can be desc ribed more 

realis tically than we do later in the simulations, see for instance the work bv lNott et al 



1999I) . they will determine, for a given value of the coefficient of restitution and in the 
hydrodynamic bulk (i.e., the region where our four hypotheses (i)-(iv) hold), which type 
of profile among those in ()5.17p - (|5.40p the system will show. 

An illustration of the phase diagram in the a-ST* plane at a given value of a < 1 is 
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Figure 3. Phase diagram illustrating the classification of Couette-Fourier fiows. This 
particular case corresponds to a = 0.9 and d = 3, as obtained from Grad's solution. 

presented in figure [S] In fact, the LTu and LTy curves have been obtained from Grad's 
solution of the Boltzmann equation (sec § U) for a = 0.9. It is apparent that the LTy 
class cannot be attained if 5T* is larger than 2y'|7(a, 0)| (~ 0.26 in the case displayed 
in figure [3]) or a is larger than aLTu(a) (— 0.36 in the case displayed in figure [3]). As the 
coefficient of restitution increases, both \j{a, 0)| and oltu decrease, so that the CTu/XTy 
and CTy regions shrink. Of course, in the elastic case only the region XTu persists. All 
these features are clearly seen in the full phase diagram depicted in figure [2] 

An interesting remark in the case of symmetric walls, i.e., ST* = 0, is the impossibility 
of having a temperature profile that is concave in the variables s or Ux but convex in 
the variable y (CTu/XTy region). As figure [3] shows, if ST* = and both plates are at 
rest (a = 0), T{y) is concave. As shearing is introduced and increased, the concavities 
of T{y) and T{ux) decrease until the value a = Oltu is reached, where the temperature 
is uniform and Ux{y) is linear (USF). Further increase of the shearing produces convex 
profiles T{y) and T{ux)- Thus, the existence of the 'hybrid' CTu/XTy region requires 
asymmetric walls {ST* ^ 0). 



6. Comparison with computer simulations 

6.1. Simulation details 

In this section we present the results obtained from DSMC and MD simulations for hard 
spheres {d = 3) and compare them with the analytical results derived from Grad's theory. 

The simulation methods that we used for DSMC and MD simulat ions are similar to those 

in our previous works and h ave been explained in detail elsewhere (^Lobkovskv. Vega Rcvcs & Urbach 
2009tlVega Reves fc Urbach 2009; Vega Rcvcs, Garzo fc San tos 2011a,: iVega Rcves, San tos fc Garzo 
201161 ) . We will briefly recall that DSMC yields an exact numerical solution of the cor- 
responding kinetic equation (inelastic Boltzmann equation in this case), whereas MD 
yields a solution of the equations of motion of the particles. Therefore, the main dif- 
ference between results from both methods is that MD simulations lack the bias of the 
inherent statistical approximation of the Boltzmann equation, where velocity correla- 
tions between particles which a re about to collide are not considered. As in our previous 
work (|Vega Reves et aZ.ll201lQl ). the global solid volume fraction in the MD simulations 
has been taken equal to 7 x 10"'^ (dilute limit), using N ^ 10^-10^ particles. In DSMC 
simulations we take a similar number of particles, N — 2 x 10^. The boundary condi- 
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System 


y'2T^/m 


T{-h/2) 


T[h/2) 


u^(-h/2)-U- 
^2T_ /m 


U+-u^{h/2) 


n(-h/2) 
n 


n{h/2) 
— 




A 


5.5 


0.9706 


7.1799 


0.1892 


0.6080 


2.1357 


0.2939 




B 


10.6 


1.2953 


8.9035 


0.2634 


0.8651 


2.8482 


0.4193 




C 


11.3 


1.3397 


9.2022 


0.2715 


0.8457 


3.0905 


0.4610 




D 


11.85 


1.3741 


9.3722 


0.2727 


0.8584 


3.1788 


0.4897 




E 


14.0 


1.5154 


10.2501 


0.2861 


0.8934 


3.5821 


0.5625 




F 


17.0 


1.7316 


10.9953 


0.3062 


0.9302 


4.1538 


0.6889 



Table 2 . Values of the wall velocity difference and of the hydrodynamic fields near the walls 
for six representative systems. In all the cases a = 0.9, h = 15(-\/27i"7^o"^)~^ and T^/T^ = 10. 



System 


K 


ST* 


a 


7 


$ 


Class 


A 


0.994 


0.1589 


0.2491 


-0.0110 


0.0218 


CTy 


B 


1.044 


0.1064 


0.3597 


-0.0022 


-0.00004 


LTy 


C 


1.038 


0.1008 


0.3697 


-0.0013 


-0.0049 


CTu/XTy 


D 


1.047 


0.0972 


0.3753 


-0.0006 


-0.0087 


LTu 


E 


1.062 


0.0854 


0.3994 


0.0017 


-0.0260 


XTu 


F 


1.073 


0.0683 


0.4251 


0.0044 


-0.0558 


XTu 



Table 3. Values of the parameters K [equation (|5.ip ]. 5T* [equation (|5.33p ]. a [equation (|3.9p ]. 
7 [equation (|4.19|l ] and $ [equation (|5.8|l ] for the systems described in tabled The right-most 
column shows the class each system belongs to. 



tions used here are analogous in both methods. When a particle collides with a wall, its 
velocity is updated following the rule v —>■ v' + C/±e^. The first contribution {v') of the 
new particle velocity is due to thermal boundary condition, while the second contribu- 
tion {U±ex) is due to wall motion. The horizontal components of v' are randomly drawn 
from a Maxwellian distribution (at a temperature I±), whereas the normal component 
vj, is sampled from a Raylc igh probability distribution: -Pd^yl) = {m\Vy\/T±)e~™''"y '"^"^^ 
( Alexander fc Garcialll997[) . 

At a given value of a, we consider a common wall distance h = \h(\/'lTfna'^^ 
n is the average density, and 8 different series of simulations with T+/T^ = 2.5, 5.0, 
7.5, . . . , 20.0. For each value of the wall temperature ratio, a number of wall velocity 
differences {U+ - [/_)/^2T_/m « 2-20 is taken. 

Once the steady state is reached, the local values of p(y), uJy), T (y) and v{y) oc 
p(y)[r(y)]~^/^ are coarse-grained into 25 layers (IVega Reves et a^.ll20116h . The local shear 
rate a is obtained from equation p.9p . Next, the local curvature parameters 7 and $ are 



15(-\/27r?icr^) ^, where 



obtained from equations (|5.6p and (|5.10p , respectively. In order to evaluate the derivatives 
dux/dy, dT/dux and d'^T/dul., the profiles Ux{y) and T{ux) are fitted to polynomials 
(typically of fifth degree) . 



6.2. Hydrodynamic profiles 

Similarly to previous works, we have observed in all simulation runs that p, a, 7 and $ 
practically remain constant in the central layers of the system. Thus, in the subsequent 
analysis the local values of p, a, 7 and $ are replaced by global values obtained by a 
spatial average in the bulk domain. 

The five classes of flows summarized in tablc[T]and figure[3]are found in the simulations. 
The USF state with thermal walls, which requires 6T* = 0, was analyzed elsewhere 
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Figure 4. (a) Profile T(y), (6) parametric plot T{ux), (c) profile p{y) = (nr/n)p{y) and (d) 
parametric plot qy{T), as obtained from DSMC simulations for the systems A (A), B (+), D 
(x), E (a) and F (■) described in table[21 Lines in (a) and (6) represent the theoretical profiles. 
Additionally, we present T{y) and T(ux) plots (O) as obtained from MD simulations for state 
E. The quantities are scaled with respect to the reference units described in the text. 





Figure 5. (a) Profile T{y) and (6) parametric plot T{ux), as obtained from DSMC simulations 
for the system C described in table[2] Lines represent the theoretical profiles. The quantities are 
scaled with respect to the reference units described in the text. 



( Vega Reyes et aZ.ll2010l l2011a ) and is not considered here. As an illustration, let us 
consider the six representative systems described in table [H We observe that, at fixed 
values h = 15{\/2TTna^)~^ and T^/T- = 10, the fluid temperatures near the walls do 
not coincide with the imposed wall values (temperature jumps). As we increase shearing, 
the differences T(±/i/2) — T± increase, changing from negative to positive values (see 
three first columns in table [2]). As for the velocity slips (jLvm 19961), i.e., the differences 
Ua:(±/i/2) — U±, they also tend to increas e (with one exception) with increasing shearing. 
In w hat follows, as in former works ( Vega Reves fc UrbachI l2009t IVega Reves et al 



2011q[ ). we take the quantities near the cold wall as reference units. Thus, m, Tr 
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T{—h/2) and Tr = l/i'{—h/2) define the units of mass, energy and time, respectively. 
Therefore, distances are measured in units of the nominal mean free path Tr^/T^jm = 
5c,,/(16Y^nrCr^); where n^ = n{—h/2). Moreover, the density is scaled with respect to 
Ur- The steady-state hydrodynamic profiles for the systems in table [2] are shown in figures 
|4]and[5j Since the profiles in system C are very close to those of systems B and D, system 
C is absent in figure |4] and its temperature profiles are shown separately in figure [3 It is 
quite apparent that the pressure is practically uniform in all the cases, thus confirming 
the hypothesis (i) made in § [5l Notice also that, even though in the simulations the size 
is fixed at h = 15(A/27rntT^)^^, the dimensionless size of systems A-E in the units of our 
choice varies since n,./n is different in each case. Moreover, in our reduced units p{y) « 1 
at all places and systems and so, for better visualization, in figure |4l^c) we choose to 
plot p{y) instead. The (bulk) temperature profile T{y) is concave for system A, linear 
for system B and convex for systems C-F. Regarding the profile T{ux), it is concave 
for systems A-C, linear for system D and convex for systems E and F. The parametric 
dependence of g^ versus T is linear (in the bulk region) in all the cases, in agreement 
with equation (j5.12p . being an increasing function for systems A-C, constant for system 
D and decreasing for systems E and F. 

The values of the quantities K, ST*, a, 7 and $ obtained from the hydrodynamic 
profiles of systems A-F are displayed in table [31 Notice in this table that the measured 
values of <& and 7 correctly predict in all cases the observed curvatures of T(y) and T{ux), 
respectively. Moreover, we have obtained a very close approach to LTy and LTu states 
in systems B and D (for which $ = —0.00004 and 7 = —0.0006. respectively). 

We introduced the simulation values of K , a, 7 and $ into the, according to our 
description, corresponding theoretical profiles for T{y) and T{ux), by using the pertinent 
(depending on the signs of 7 and $) expressions given in ij 15.31 It is worth remarking that 
the theoretical profiles T{y) do not depend on the separate values of K , a, 7 and $ but 
only on the two combinations Tq and Iq [cf. equations ()5.15p ]: as for the theoretical profiles 
Tiux), they depend on the same parameter Tq as before plus the combination aw. The 
resulting profiles are included in figurcs|4lja),|4|6) and El where the integration constants 
2/0 and uq are determined as to reproduce T and u^^ at y = 0. As we can observe, the 
agreement between the theoretical curves from our generalized hydrodynamic description 
(see § 15. 3p and simulation data is excellent, the deviations typically being restricted to 
1-2 layers near the cold wall and 2-4 layers near the hot wall. Those small deviations 
can be due to boundary-layer effects and/or to residual limitations of the hydrodynamic 
description exposed in § [51 In any case, it is worth remarking that the local mean free path 
(inversely proportional to the local density) is larger near the hot wall (where deviations 
present a longer range) than near the cold wall. As a matter of fact, in the employed 
reference units, the mean free path is ^ 1 near the cold wall and ^ n{—h/2)/n{h/2) = 6- 
7 near the hot wall. It is also interesting to note that the lack of agreement near the 
boundaries seems to become less important as the shearing increases (i.e., from system 
A to system F). 

The simulation data plotted in figures IH and [H have been obtained from the DSMC 
method but they perfectly agree with those obtained from MD. As an example, we 
compare the results obtained from both simulation methods in one of the curves of 
figures mj a) and 111^6). 

6.3. Transport coefficients 

Once we have checked that the steady base states discussed in § [5l are supported by the 
simulations, we now proceed to present the simulation results for the transport coefficients 
and compare them with Grad's theoretical predictions. 
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Figure 6. Thermal curvature parameter 7 as a function of shear rate squared a^ for two values 
of the coefficient of restitution: (o) a = 0.9 and (6) a = 0.7. Lines represent results from Grad's 
analytical solution (solid lines) and from the NS prediction (dashed lines), while symbols stand 
for DSMC (□) and MD (■) simulations. 



As a general trend, we have observed a relatively good semi-quantitative agreement 
between simulation and Grad's theory for all relevant quantities, except for the reduced 
thermal conductivity A* and for the reduced viscosity 77* at low a. In figurelHlwe plot the 
results for the thermal curvature parameter 7 for two different values of the coefficient 
of restitution; a = 0.9 and 0.7. We detect, both in simulations and theory, the aforemen- 
tioned transition from 7 < for low shear rate s to 7 > for higher shear rates. This 
transition is also predicted by the NS solution (IVe ga Reves &: UrbachI 120091 ). in which 
case 7 is a linear function of a^ [see equation p.l2p ]. As we see, the true parameter 7 
has a more complex dependency on a. It is apparent that Gra d's theory predicts well 
the value a = oltu where 7 = 0, as already shown elsewhere ( Vega Reves et aLll2010l . 



2011q| ). It is also noteworthy that, in the region 7 > 0, Grad's theory does a better job 
for a = 0.7 than for a ~ 0.9. It might seem surprising that both NS and Grad's pre- 
dictions for 7 show significant discrepancies with simulation data in the region of small 
shear rates, especially for a = 0.7. The explanation lies in the fact that, apart from a 
and ST*, 7 is an additional measure of the strength of the gradients, which in the limit 
a ^ is governed by a and thus cannot be done arbitrarily small for finite inelasticity. 

As discussed in § 15.31 for a given value of a, it is possible to find pairs ((5T*, a) such that 
the temperature profile T{y) is linear (LTy states). It is also possible to find a value of 
a, independent of ST* , where the temperature proffie T{ux) is linear (LTu states). These 
two loci split the plane ST* vs a into the three regions sketched in figure [3] We represent 
in figure [7] the phase diagram, as obtained from our simulations, for (a) a = 0.9 and (6) 
a = 0.7. For comparison, the curves predicted by Grad's solution are also included. As we 
see, the agreement between theory and simulation is qualitatively good for both values 
of a. As a complement, figure [8] shows the threshold value a^rp versus the coefficient of 
restitution for ST* = 0.015. We observe that the LTy is not possible for this value of the 
slope ST* if a ^ 0.967. 

In figures l9l and [TOl we plot the shear-rate dependence of the reduced shear viscosity 77* 
and of the normalized diagonal components of the stress tensor 9^ and 0y, respectively. 
It is quite apparent that, except for the shear viscosity in the range of low shear rates, 
the agreement between Grad's analytical solution and DSMC and MD simulations is 
quite good (somewhat better for a = 0.9). The agreement is specially good around the 
LTu states (i.e., a^ y 0.15 and a^ ?^ 0.55 for a = 0.9 and a = 0.7, respectively), as 
previously reported (jVega Reves et aZ.ll2010ll2011a[ ). Figure |9] shows that the non-linear 
shear viscosity decreases with increasing shear rate ('shear thinning' effect). In what 
concerns the reduced directional temperatures, figure [10] shows that 9x (Oy) increases 
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Figure 7. Phase diagram in the plane 5T* vs a (see figures [2] and [3} for two values of a: (a) 
a = 0.9 and (6) a = 0.7. Lines stand for the analytical solution from Grad's method, while 
open and solid symbols stand for DSMC and MD simulations, respectively. The LTy curve is 
represented by solid lines (theory) and squares (simulation), while the LTu line is represented 
by dashed lines (theory) and triangles (simulation). 




Figure 8. Threshold value flLTy for which the linear T{y) occurs if ST* = 0.015, as a function 
of the coefficient of restitution. Line stands for Grad's method solution, while open and solid 
symbols stand for DSMC and MD simulations, respectively. 



(decreases) with increasing shearing. It is interesting to note that 9x < 9y for very 
small shear rates, until both quantities cross at a certain value of a. This phenomenon 
is qualitatively captured by Grad's solution. Comparison between figures [Hlja) andj^Jfe) 
shows that, as the inelasticity decreases, the region of shear rates corresponding to 7 < 0, 
and hence the region with worse Grad's predictions, shrinks . In fact, in the purely elastic 
case {a = 1) the Grad expression for rf is rather accurate (JGarzo fc Santosll2003f ). 

Finally, in figure [TT] we plot the results for the two heat flux transport coefficients 
(thermal conductivity A* and cross coefficient </>*). As already explained, there is in gen- 
eral a (non-Newtonian) horizontal component of the heat flux, from which the cross 
thermal conductivity coefficient cj)* results. Perhaps surprisingly, we find that the agree- 
ment between Grad's theory and simulations is better for the cross coefficient 4>* than 
for the thermal conductivity A*. Moreover, while Grad's theory predicts that A* weakly 
increases with a (a = 0.9) or exhibits a non-monotonic behavior (a = 0.7), simulations 
yield a decreasing A* vs a. On the contrary, the agreement for the cross coefficient is 
qualitatively good, since ip* vs a is increasing both for Grad's theory and simulation. 
This agreement is very good in the region of low shear rates up to the threshold value 
for LTu states (as expected), whereas for higher shear rates the theory and simulation 
results tend to separate. 

A final comment regarding the comparison between simulation and Grad's theory is in 
order. According to equation (|4.18p . (* ex k'r]*a'^ - X*j, where k = {d-l)/d{d + 2). Since 
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Figure 9. Generalized viscosity 77* as a function of a^ for (a) a. = 0.9 and (6) a = 0.7. Lines 
stand for Grad's method solution, while open and solid symbols stand for DSMC and MD 
simulations, respectively. 
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Figure 10. Reduced normal stress components 9x (dashed lines and triangles) and 6y (dot- 
ted lines and squares) as functions of a^ for (a) a — 0.9 and (6) a = 0.7. Lines stand for 
Grad's method solution, while open and solid symbols stand for DSMC and MD simulations, 
respectively. 
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Figure 11. Generalized thermal conductivity A* (solid lines and triangles) and heat flux cross 
coefficient 0* (dashed lines and squares) as functions of a^ for (a) a = 0.9 and (b) a — 0.7. 
Lines stand for Grad's method solution, while open and solid symbols stand for DSMC and MD 
simulations, respectively. 



the reduced cooling rate Q* is satisfactorily captured by Grad's method [(see equation 
p.7p ] , we conclude that the deviations of if , A* and 7 from the simulation data are not 
entirely independent and are somewhat constrained by the combination j^ri*a^ — A* 7 
(note that k = jf for d = 3). In fact, figures |6l l9l and [Til show that, in the region with 
7 < 0, I7I and A* are underestimated by Grad's solution, while -q* is overestimated. In 
the region of 7 > 0, however, ry* is quite accurate, so that the underestimation of 7 is 
compensated by an overestimation of A*. It is interesting to remark that the accuracy 
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Features 


Level of description 








NS Grad Gener. non-Newton. 


Simulation 




(i) p = const 


Derived Assumed Assumed 


Observed 




(ii) a — const 


Derived Derived Assumed 


Observed 




(iii) P,, ^ F{dyT) 


Construction Derived Assumed 


Observed 




(iv) Qy ocdyT 


Construction Derived Assumed 


Observed 




Pxx 7^ Pyy J^ Pzz 


No Yes Yes 


Yes 




q.^O 


No Yes Yes 


Yes 




Transport coefficients 


Explicit Explicit Unspecified 


Measured 



Table 4. Hypotheses (i)-(iv) and main features of the plane Couette-Fourier fiow according to 
the level of description; NS (§[31, Grad's 13-moment method (§[4|, generalized non-Newtonian 
hydrodynamics (§[5]) and simulation (§|6]). 



of Grad's quantitative predictions is highly correlated with the magnitude of the ther- 
mal curvature parameter 7, i.e., the smaller I7I the better the general performance of 
Grad's solution. In fact, the agreement between theory and simula t ion is quite g ood in 
the LTu state (7 = 0), as previously shown by IVega Reves etal\ (|2010l |2011q() . This 
confirms the role played by 7 as an intrinsic measure of the strength of the gradients 
(IVeea Reves fc Urbach|[20oi . 



7. Conclusions 

7.1. Summary 

We have studied in this paper the laminar flows in a low density granular gas conflned 
between two infinite parallel walls, which, in general, are at different temperatures. Ad- 
ditionally, the granular gas can be sheared if there is relative motion between both walls. 
We have described a general classification of steady granular Couette-Fourier flows that 
occur in this system, at constant pressure, for arbitrarily large velocity and temperature 
gradients. We have shown that, due to symmetries in the system, the steady-state equa- 
tions for the flow velocity and temperature are quite simple, even in the non-Newtonian 
regime, and have a straightforward analytical solution. Moreover, the type of solutions 
for the hydrodynamic profiles turn out to be dependent on just two constant parameters: 
the thermal curvature coefficients 7 and $. The former is proportional to the second 
derivative of T in a spatial variable scaled with collision frequency, while $ is related to 
the second derivative in the natural spatial variable. Depending on the different possible 
combinations of signs of these two parameters, the corresponding steady profiles can be 
grouped into five different classes of flows, each one having peculiar properties (see table 

The main conclusions of this work are that the assumptions made on the form of the 
hydrodynamic profiles [see equations p.9|) and (j4.13p - (|4.16|) ]. as well as the associated 
classification of fiows, have been validated by three independent routes. From a theoretical 
perspective, we have obtained an exact solution of the set of moment equations derived 
from Grad's method applied to the inelastic Boltzmann equation. Next, we have simu- 
lated the Couette-Fourier flows by using the DSMC method (which numerically solves 
the Boltzmann equation) and MD simulations (which numerically solve the equations of 
motion of the system of inelastic hard spheres). 

This triple validation extends in a non-trivial way some of the qualitative features of 
the NS description to the realm of non-Newtonian hydrodynamics. This is summarized 
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in table IH As shown in § [3l the NS constitutive equations, complemented by the mo- 
mentum and energy balance equations in the steady Couette-Fourier geometry, imply 
the fulfillment of points (i)-(iv) without further assumptions. On the other hand, they 
do not account for normal stress differences or a heat flux component parallel to the 
flow. This is remedied by Grad's moment method, in which case only hypothesis (i) on 
the constancy of pressure is needed. A more general non-Newtonian treatment makes 
use of the four assumptions on the same footing, thus allowing us to accommodate for 
any specific form of the generalized transport coefficients. Finally, the simulation results 
are seen to support the validity of those assumptions, providing as well the dependence 
of the main quantities on both the shear rate and the coefficient of normal restitution. 
However, it must be kept in mind that, while simulations are essentially consistent to 
a large extent with the generalized hydrodynamic description of § El slight deviations 
due to the high intricacy of the Boltzmann equation cannot be discarded. Those small 
de viations have been reported in the case of th e pure Fourier flow for elastic hard spheres 
by iMontanero. Alaoui. Santos fc Garzd ( 1994 ) . 



While Grad's moment method supports the four assumptions (i)-(iv), as well as the 
existence of normal stress differences and a heat flux component orthogonal to the thermal 
gradient (see table S]) , we have observed that a quantitative agreement with simulations 
is generally good near the LTu state (i.e., for small values of I7I) only. As the magnitude 
of the thermal curvature parameter 7 increases, some transport coefficients (77* for 7 > 0, 
(j)* for 7 < and 9j. and 9y for both 7 < and 7 > 0) are better predicted by Grad's 
theory than other ones (ry* for 7 < 0, (/)* for 7 > and A* for both 7 < and 7 > 0). 

7.2. Discussion 

The signs of 7 and $ depend on both the physical properties of the granular gas and 
the boundary conditions. However, rather than analyzing the interaction between gas 
and wall, our work is f ocused, similarly to previous works ( Vega Reves fc UrbachI 120091 : 



Vega Reves et al\\201(t , on the bulk properties of the gas itself and we study all possible 



transitions between the different classes of flows. All class transitions have been generated 
by using the usual hard wal l boundary conditions, both in D SMC and MD simulations 
(see for instance the work bv lGalvin. Hrenva fc Wildmanll2007l where the same boundary 



conditions are used for simulation of thermal walls) . The phase diagram obtained from 
simulations is completely analogous to the theoretical one, depicted in flgures ^ and [31 
as shown in flgurejTl We have checked in the simulations that, as in figure |21 only two 
of the five possible flow classes (see table |T]) define surfaces in the three-parameter space 
{a, 6T*,a}. They divide this space into three regions that define three other entire classes 
of granular flows. Thus, we have taken these surfaces as a reference for our analysis of 
flow class transitions. One of the surfaces is the LTu flow class (7 = 0), characterized 
by linear temperature vs flow velocity profiles, and already studied in former works 
( Vega Reves et aZ.ll2010ll2011a ). The other surface is the LTy class ($ — 0), characterized 



by linear temperature vs vertical coordinate profiles. The LTy surface is always below 
(lower shear rates) the LTu surface (figure H]), except for walls at the same temperature 
(ST* — plane), where they coincide, defining a curve that is the remaining sixth flow 
class, which can be regarded as a subclass of the LTy or LTu classes. This class (or 
subclass) is actually the well known uniform shear flow (USF), i.e., constant T and linear 
Ux{y)- Note that here the USF i s achieved with therma l walls rather than with generalized 
periodic boundary conditions ( Lees fc Edwards! 1 19721 ). Regarding the other classes, the 



first region (CTy) is below the LTy surface and is characterized by 7 < and 4> > 0. 
The second region (CTu/XTy) occupies the space between the LTy and LTu surfaces, 
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being characterized with 7 < and $ < 0. FinaUy, the third region (XTu) is above the 
LTu surface and corresponds to 7 > and $ < (see figure [2]) . 

One important difference between LTy and LTu classes is that, while LTu flows are 
possible for arbitrarily large 6T*, the LTy flows are restricted to values of ST* smaller 
than a threshold value (5T£rp (a, a), which has an upper bound at a = (see figures[3]and 
[7]). The agreement between theory and simulation in this aspect is qualitatively good. 
In particular, we have checked that a too large ST* in the simulations results in a direct 
LTu transition without passing through an LTy transition, when increasing shear rate 
from a = 0. For instance, for a ~ 0.9, and following results in figure [7Ka), a value of 
ST* = 0.3 suffices to suppress the LTy transition. Thus, in this case we would already 
start from $ < at a = 0, never entering the class of flows with concave T{y). 

We have not detected so far instabilities (departures from laminar flows) in the simu- 
lations. This is reasonable since the flows that we have analyzed are either below or not 
far above the LTu surface, and thus the y occur at low Reynolds number Re (LTu flows 
typically have Re ^ 100, see the work bv lVega Reves et al\\2011m . In order to see higher 



Re we would need to separate much further above the LTu surface, at extremely large 
shear rates, or apply larger ST*. 

In conclusion, we have described in detail, by means of theoretical and computational 
studies, all possible classes of base laminar flows for a low density granular gas in a 
Couette-Fouricr flow geometry. Those classes differ in the curvature of the T{y) and 
T{ux) profiles but otherwise they can be described within a common framework char- 
acterized by a heat flux proportional to the thermal gradient and uniform stress tensor 
and reduced shear rate. This unified setting encompasses known and new states, from 
the Fourier flow of ordinary gases to the uniform shear flow of granular gases, from the 
symmetric Couette flow of ordinary gases to Fourier-like flows of granular gases with con- 
stant thermal gradient and from states with a magnitude of the heat flux \q\ increasing 
with temperature to states with a decreasing, a constant or even a zero |<7|. 



7.3. Outlook 

The flow classes described in this work might be useful for future works in a variety of 
problems in granular dynamics , such as the study of a. granul ar impurity under Cou- 
ette flow (JGarzo fc Vega Revesll2010t IVega Reves et aUuQWm . This implies that the 
same set of flow classes should exist for the granular impurity; LTu and LTy classes for 
instance. This may have implications to segregation conditions for a granul ar impurity 
(| Jenkins fc Yoonll2002HGarz6 fc Vega RevedbooJICarzo fc Vega ReveslboiOl ). Moreover, 
a comp lete determination of the steady base states is convenient for studies of insta- 
bilities (IHopkins fc Louge 1991 : Wang. Jackson fc SundaresanI 1996 : Alam fc Nott 1998 



Nott. Alam. Agrawal. Jackson fc SundaresanllQQQtlKhain fc Meersonl2003tlAlam. Shukla fc Ludine 

20081 ). Furthermore, analogous temperature curvature properties are observed for the 

same geometry in moderately dense granular gases, except that for higher dens i ties a 

region with temper ature curvature inflections grows from the boundaries ( LunI Il996t 

Alam fc Notti 119981 ). Thus, we expect some of the conclusions of the present analysis 

to be useful for instability in quite generic problems of granular flow. We are currently 

working on extensions of this work in granular segregation and flow instability. 
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Th e expressions for the NS transport coefficients are (jBrev et a/.lll998l : iBrev fc Cubero 

1 



Here, 






AiNs(") 



"/3i(a) + iC(a)' 

1 
/32(«)-^C(a)' 



/3i(a) 



/32(a) 



/32(a) - 5^C(a)J [/32(a) - 2(i^C(a) 
1 + a 



2 

1 + a 



-y<'-' 



3d + 8, 

IH 1-a 



(Al) 
(A 2) 

(A3) 

(A 4) 
(A 5) 



and the reduced coohng rate C*(a) is given by equation (j3.7p . In equations (j3.7p and 
(|A1[) - (|A5|) . terms associated with the deviation of the homogeneous coohng state dis- 
tribution from a Maxwelhan have been neglected (jGarzo. Santos fc Montanerdl2007l ). 



Appendix B. Explicit expressions in Grad's approximation 

Taking into account in equations (|4.7p - (|4.1ip the form of the fluxes given by equations 
(|4.13p - (|4.16p . one gets, after some algebra, 

2d 



a[9y-{l3^+CW] 



d-1 



2d 



70* = 0, 



(/3i + C)Ox - 2rfa' + ^— 7A* = /3: 



d-1 

6d 
d~l 



Wi + nOy + :r^iX* = Pi, 



{d + 4)a 



V 



d-1 
d + 4:„ d+2 



-X* 



/32A* 



- (d + 2)/?2r = 0, 



(Bl) 

(B2) 
(B3) 
(B4) 

(B5) 



2 ^ 2 "'"■■ ' d-1" 

The algebraic equations (jB 1[) - (|B5[) allow one to express 77*, A*, Ox, dy and </>* in terms 
of a, a and 7 as 

jf = A-i {2d^{d + 4)/3ia2 -{d- lf{d + 2)2/3i/3| + 2d [d{d + 4) {{d + 4)/3i - 2^^) 

-6{d- l){d + 2)13-2] 1}, (B6) 

A* = A-\d-l) { [2^1 -{d + 4)/3i] [{d - l){d + 2)^i/32 + 2d{d + 4)7] - 2d(d + 4)/3ia2} , 

(B7) 
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9^ = {AP^)-'^ll3i{2a^ +'p\) [2d?{d + 4)a^ -{d-lf{d + 2fl3l] +2d[2d{d + A)a^ 
X {{d + 2)/3i - 2^1) -{d-l){d + 2)h (12a' + ^1 (2^i + 3(rf + 4)/3i))] 7 
-8d2(d + 4)[^i + (rf + 4)/3i]72}, (B8) 

6iy = A-i{2rf2(d + 4);5i/3ia' - [(d - \){d + 2)/3i/32 + 12^7] 

X [(d - l)(d + 2pi/32 + 2d{d + 4)7] }, (B 9) 

i>* = A-i(d - \){d + 4)a {d^i [2;5i - (d + 4)/3i] - (d - l)((i + 2)l\^2 - 12^7} , (B 10) 
where /S^ = /^i + C* and 

A = 2d^{d + 4)(^i - 67)0' - (d - l)2(d + 2)2^i/3| - M{d - \){d + 2)(d + 4)^i/327 
-12^2(^ + 4)272. (Ell) 

Finally, substitution of rf and A* into equation (j4.18p yields a quadratic equation for 
7. Its physical solution gives 7 as a function of the shear rate a and the coefficient of 
restitution a. 

Setting 7 = in equations (|4.18p . (|B6P and ()B7p . we get the prediction for the LTu 
threshold shear rate in Grad's approximation. The result is 

aLTu(«) - J^^i. (B12) 

The expressions for the LTu transport coefficients rf ^ A*, Qx-, Oy and (/)* are obtained by 
making a = a^,tu and 7 = in equations (IB 6[) - (|B 11[) . The explicit expressions have been 
given elsewhere (jVega Reves et aLll2011af ). 

In the absence of shearing (a — ?► 0), equations (|4.18p and (|B6|) - (|B lip yield 
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+dC [3(d + 2)(d - l)/32 + d(d + 4)C*] }/(d - l)(d + 2)2/32^', (B 14) 
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(.-l)(. + 2)3,2^?(c^+2)2,, + 2(.2+3.-2)C-{^(^^^)^^^+^^[(^^^)^ 
X (d - l)/32 + 2d\2d + 7)C] + dC [3(d + 2)(d - l)/32 + (3(i2 - lOd + 4)C*] }. 

(B18) 
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In the elastic case {(* -> 0, /3i -> 1, (32 -^ 1), one has 9x -^ 1, 9y ^ 1, A* -> 1, 77* -> 1, 
7 -)► and <f>*/a -J> (2d - l)(d + 4)/((i - l)(d + 2), which corresponds to the Fourier flow 
of conventional gases. 
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